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Using molecular dynamics simulations we study supercritical fluids near the gas-liquid critical 
point under heat flow in two dimensions. We calculate the steady-state temperature and density 
profiles. The resultant thermal conductivity exhibits critical singularity in agreement with the 
mode-coupling theory in two dimensions. We also calculate distributions of the momentum and 
heat fluxes at fixed density. They indicate that liquid-like (entropy-poor) clusters move toward 
the warmer boundary and gas-like (entropy-rich) regions move toward the cooler boundary in a 
temperature gradient. This counterflow results in critical enhancement of the thermal conductivity. 
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I. INTRODUCTION 

In one-component fluids, the density and energy fluc- 
tuations are known to become long-ranged and long-lived 
as the temperature T and the density n approach the crit- 
ical values T c and n c 0, ■ The critical singularities are 
characterized by the correlation length £, which grows 
as £q(T/T c — l)~ v on the critical isochore with £ being 
a microscopic length and v being the critical exponent. 
On one hand, the isothermal compressibility Kt, the iso- 
baric thermal expansion coefficient a p , and the isobaric 
specific heat C p grow as £ 2 with fj being the small 
Fisher critical exponent. On the other hand, the thermal 
diffusivity Dt behaves as kBT/r]^ d ~ 2 and the life time 
of the critical fluctuations grows as T£ = £ 2 / Dt ~ 
where d is the space dimensionality and the weak singu- 
larity of the shear viscosity r\ is neglected. As a result, 
the thermal conductivity A = DtC p grows as £ 4 ~ !i ~'? 
The critical behavior of A and 
by the mode-coupling theory I 
renormalization group theory [5j. 

However, the calculations in these dynamical theories 
are performed in the space of the wave vector of the fluc- 
tuations and are rather formal. The real space picture 
of the enhanced heat transport in a small temperature 
gradient dT/dz is as follows 2\. The critical fluctua- 
tions with relatively higher (lower) densities should be 
convected in the direction (reverse direction) of the tem- 
perature gradient. The typical velocity of the clusters 
with lengths of order £ is given by 



has been well described 
, and by the dynamic 



move only over the distance 

An 

n 
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where An = a p LdT / dz is the difference between the den- 
sities at the two ends of the cell. This distance is very 
short for L much longer than £. Hence it should be diffi- 
cult to unambiguously observe the the cluster motion in 
heat flow experimentally and even numerically. 

As numerical work of heat conduction in fluids 0, 0] , 
the thermal conductivity has been calculated using equi- 
librium MD simulations 0, Q on the basis of the Green- 
Kubo formula | 2i | q or using nonequilibrium MD sim- 
ulations [lOl fill Il2| . In particular, developing a sim- 
ple method, Ohara performed nonequilibrium MD sim- 
ulations for Lenard- Jones (LJ) fluids 13] and for liquid 
water 0. All these previous papers treated fluids far 
from the critical point. In this paper we will use Ohara's 
method to realize heat-conducting states in the one-phase 
region near the critical point. 

This paper is organized as follows. In Section 2 we will 
present numerical results on equilibrium critical behavior 
in supercritical LJ fluids. In Section 3 we will show nu- 
merical results on near-critical heat conduction together 
with theoretical interpretations. In particular, we will 
confirm the cluster convection mechanism by introduc- 
ing steady-state distributions of the momentum and heat 
fluxes at fixed density. In Appendix B we will summarize 
the mode-coupling theory for the thermal conductivity. 
In Appendix C we will examine the linear response to 
heat flow and justify Eq.l. 
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in the linear response with /3 being the critical exponent. 
The entropy of the liquid- like regions is smaller than that 
of the gas-like regions by n^ d Ss ~ where the 

entropy fluctuation Ss (as well as the denisity fluctuation 
Sn) has sizes typically of order (~ (T/T c - if 

on the critical isochore). The thermal average of the 
convective heat flux nTSsv^ thus gives rise to the critical 
heat conduction. Within their life times the clusters can 



II. MODEL AND EQUILIBRIUM RESULTS 

We used a two-dimensional (2D) LJ fluid composed of 
N identical particles. The pair potential as a function of 
the distance r between two particles is given by 



(r) = 4e 



r , 



a-- 1 - 
r 



C (r<r c ), 



(3) 



and 4>(r) = for r > r c . The constant C is chosen such 
that 4>(r c ) = 0. The cut-off length r c was set equal to 3<r. 
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TABLE I: Numerical estimations for the two-dimensional 
LJ fluids. The estimated critical temperature T c and density 
n c are given in the first and second columns. The particle 
number N and the cut-off radius r cut used are given in the 
third and fourth columns. 
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2.5 


[17] 


0.498 


0.36 


8000 
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0.47 


0.35 


4096 
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[19] 


0.47 


0.37 


5000 


3.0 


Present work 



"The Gibbs ensemble method was used with r cu t 
is the cell size. 



L/2, where L 



The system contains N — 5000 particles. Space and time 
are measured in units of a and To = (ma 2 /e) 1 / 2 , where m 
is the particle mass. Equilibrium states of the fluid may 
be characterized by the temperature T and the average 
number density n — N/V measured in units of e/fce and 



a , respectively. The pressure is measured in units of 



a €. We used the leapfrog algorithm to integrate the 
Newton differential equations with a time step of 0.01 
and the cell-index method of cutting off the interaction 
potential. The details of these numerical methods are 
described in the literature Q. 

The phase diagram of the two-dimensional L J fluid has 
been studied by several groups using the conventional 
Monte Carlo method the Gibbs ensemble method 
EU El 113 j an d finite-size scaling analysis . Table |I] 
summarizes the critical parameters reported in the liter- 
ature together with our MD results. However, note that 
the critical parameters largely depend on the details of 
the truncation of the potential [17|. We also mention 
that Luo et al [2(j examined thermal relaxation in a two- 
dimensional supercritical LJ fluid. 

As preliminary work before nonequilibrium simula- 
tions, we carried out equilibrium simulations in the 
canonical (constant-iVV'T) ensemble, using the Nose- 
Hoover thermostat 0, El under the periodic bound- 
ary condition, with At = 0.01, to calculate the struc- 
ture factor S(q). We started with random initial parti- 
cle configurations at each given temperature, waited for 
t w = 5 x 10 4 , and afterward took data in a subsequent 
time interval of t w < t < 2t w . This long equilibration is 
needed because the density fluctuations relax very slowly 
near the critical point. That is, i w should be longer than 
the life time = £ 2 /Z>r of the critical fluctuations with 
sizes of the order of the correlation length £ Q. Here A 
is the thermal conductivity and C p is the isobaric heat 
capacity per unit volume, respectively. In particular, in 
two dimensions the critical exponent v is equal to 1 and 
the critical singularity of Dt is weak (as will be evident 
in Eq.19), so on the critical isochore grows as 



n ~ e ~ (t/t c 



1 



(d 



(4) 



We consider the structure factor given by 

5(g) = f dre l ^ r (n(r,t)n(0,t))/n. (5) 



Here we define the fluctuating particle number density in 
terms the particle positions as 



Kr,i) = ^£(n(*)-r). 



(6) 



i=l 



We took the angle average in the calculation of S(q). 
Fig.l shows S(q) for T = 0.65, 0.51, 0.5, 0.495, and 0.49 
at n = 0.37. We can see the power law q behavior, 



S(q) ~ q~ V \ 



(7) 



in the range £ -1 < q < 2 near the critical point, where 
the exponent value 7/4 is consistent with the well-known 
Fisher critical exponent r) = 1/4 in two dimensions. The 
peak around q ~ 6 represents the short-range pair cor- 
relation at this density. We then determined the corre- 
lation length £ by fitting the data to the extrapolated 
expression S(q) = nk B TK T /[l + (q0 2 } 7/8 for q < 1. 
The isothermal compressibility Kt = (dn/dp)? /n can 
be determined from the long wavelength limit of S(q). 
We show £ vs n in Fig. 2 and Kt vs n in Fig. 3 for var- 
ious T(> 0.49). Although not shown in Figs. 2 and 3, 
we also performed simulations at lower temperatures to 
obtain £ — 50 for T = 0.485 and £ > L (apparently) 
for T = 0.48 at n = 0.37. However, for these T, our 
simulation times are not sufficiently long compared with 
t c in Eq.4. From the peak positions in Figs. 2 and 3 we 
estimated n c = 0.37 in Table 1. This value was also ob- 
tained as the mean position of the two peaks in the one- 
body density distribution ^>(p) (defined by Eq.22 below) 
in equilibrium for 0.495 < T < 0.50. 

From the data in Figs. 2 and 3 on the critical isochore, 
Kt behaves as a function of £ as 



K T = 3.7£ 7/4 



0.80, 



(8) 



in units of a 2 /e. We then fitted £ to the scaling form 
£ = £o(T/T c — 1) _1 on the critical isochore to obtain 
T c = 0.47 and £o = 0.6. From the isothermal curves in 
the p — n plane in the range 0.495 <T< 0.50 at n = n c , 
we obtained 221 



dp 
df 



dp 
df 



— [ ] = 0.40. 



(9) 



We also consider the specific heat C p = nT(ds/dT) p per 
unit volume (s being the entropy per particle) and the 
thermal expansion coefficient a p — —(dn/dT) p /n at con- 
stant pressure. These quantities grow strongly and are 
related to Kt by 



K i 



(10) 
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near the critical point. These relations will be used in 
the next section. 

As long as £ <C L, our equilibrium results are consistent 
with the well-known results of critical phenomena 0, Q • 
If £ approaches L, the finite-size scaling analysis may be 
performed [l9|. However, such analysis is beyond the 
scope of this paper. 

III. NONEQUILIBRIUM SIMULATIONS 
A. Method 

Next we imposed a heat flux to the system using 
Ohara's method 0, As illustrated in Fig.4, the 

cell is divided into three parts, cooling, heating, and in- 
terior regions. In the cooling region — 0.5L < z < — OAL 
the average temperature of the particles was kept at Tl, 
while in the heating region OAL < z < 0.51/ it was kept 
at Tjj- The precise definition of the average temperature 
in a given region will be presented in Eq.12 below. The 
pinning of the average temperatures in the cooling and 
heating regions was realized by simple scaling of the ve- 
locities of the particles in the two regions at every time 
step. The periodic boundary condition was imposed in 
the x direction, while the walls at z = ±L/2 were as- 
sumed to interact with the particles via the LJ poten- 
tial in Eq.3 where r is the distance from the wall and 
r c = 3. The particles in the interior (— OAL < z < OAL) 
obeyed the Newtonian dynamics without artificial ther- 
mostat. The particles entering the interior from the cool- 
ing (heating) region have lower (higher) kinetic energies 
than those of the particles in the interior on the average. 
Then a steady heat conducting state is realized after a 
transient time. 

In our nonequilibrium simulations, we used a single 
density n = 0.37 nearly equal to n c . The system length 
is then L = (5000/0.37) 1 / 2 = 116. The lower boundary 
temperature Tl was changed as 0.7, 0.65, 0.6, 0.52, 0.505, 
0.5, 0.495, and 0.49. The temperature difference AT = 
Th — Tl was fixed at 0.005 in all the simulations. We 
regard the system to be in a steady state for t > t w — 
6 x 10 4 after application of AT. In the following the 
steady-state values of the physical quantities will be the 
time averages over the data during the next time interval 
t w <t <t w + < data with i data = 14 x 10 4 . 

B. Steady-state density and temperature profiles 

Fig. 5 displays a snapshot of the particle positions in 
the cell at t = 2 x I0 5 for T L = 0.50, where the sys- 
tem is nearly in a steady state. The large clusters 
formed by many particles are significantly denser near the 
cooler boundary (bottom) than near the warmer bound- 
ary (top) [23 • This is due to the diverging isobaric ther- 
mal expansion as will be shown in Eq.16 below. By com- 
paring successive snapshots (not shown here), we recog- 



nize that the clusters appear and disappear continuously 
on the time scale of T£ in Eq.4. 

To quantitatively analyze Fig. 5, we need to calculate 
the time averages of the temperature and the density. 
They are defined as follows. Dividing the interior into 
eight layers with thickness L/10, the density in the £-th 
layer is defined by ni(t) = (10/ L 2 )Ni{t) in terms of the 
particle number in the £-th layer, 

N e {t) = / dz dxn(r,t), (If) 

Jz e JO 

where zt = (£ — 5)T/10 and n(r,t) is the fluctuating 
density in Eq.6. The temperature in the £-th layer Tg(t) 
may be defined by 

m) = j^J2\vi(t)-Mt)\ 2 , (12) 

where the summation is over the particles within the £-th 
layer and vg(t) = Yliet V i(t)/Nt(t) is the average velocity 
within the £-th layer. Notice that the 0-th layer is the 
cooling region and the 9-th layer is the heating region. 
Thus we set T (t) = T L and T 9 (t) = T H in the cooling 
and heating regions, respectively. 

Fig. 6 shows the steady-state temperature and density 
profiles, 

T(z) = (T i (t)), n(z) = <»*(*)>> (13) 

which are the time averages with z — (£ — 5)L/10. Here 
Tl = 0.50 and n — 0.37 as in Fig. 5. For this case the 
deviation of T(z) from the linear profile is not large and 
there exists a temperature gradient also in the cooling 
and heating regions. We may define the penetration 
widths and dn by the extrapolations, 

T(-d L -2L/5) = T L , T(2T/5 + d H ) = T H . (14) 

When the temperature profile is nearly linear, we simply 
find c?l — c?h — L/20. When the temperature flux is 
not too large, the effective cell length of heat conduction 
becomes 

4L/5 + rf L + c?h = 9T/10, (15) 

which is the distance between the middle points of the 
cooling and heating region as in the case of Ohara's sim- 
ulation [iH . 

On the other hand, in Fig. 6 the density deviation is 
much more enhanced than that of the temperature. We 
expect that if the deviation 8T(z) = T(z) — T(0) (mea- 
sured from the center y = 0) is not too large, the average 
density deviation Sn(z) = n(z) — n(0) should be given by 

5n(z) = ~na p 8T(z), (16) 

where a p is the isobaric thermal expansion coefficient. 
Here, at the center, we find T(0) = 0.5025 and n(0) = 
0.37, so Eq.10 indicates Kt = 159 and na p — 23.5 
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at the center. In Fig. 6 the solid line has a slope of 
AT/(9L/10) = 0.0056/L, while the dotted line a slope of 
—23.5 x 0.0056/L. These two lines can fairly fit the tem- 
perature and density data, though there are considerable 
deviations close to the boundaries (obviously because the 
boundary regions are considerably off-critical). Notice 
that we assume homogeneity of the pressure in Eq.16. 
To check this, we calculated the steady-state time aver- 
ages of the trace of the stress tensor (integrated in each 
layer) and found no appreciable heterogeneity in these 
values. 



C. Thermal conductivity 

In our simulations, the steady-state thermal conduc- 
tivity A in the interior was calculated from 



X = Q 



0.9L 
AT ' 



(17) 



where 0.9L is the effective cell length in Eq.15. The Q is 
the steady-state heat flux written as 



Q = -(J Qz (t))/0.8L 2 



(18) 



where Jq it) is the integral of the heat flux density 
within the interior with area 0.8L 2 and its microscopic ex- 
pression will be given in Appendix A. In our small system 
the fluctuations of the heat flux turned out to be large, 
so we performed 10 independent runs and calculated the 
mean values of the corresponding 10 time averages. In 
Fig. 7 the thermal conductivity data are shown as a func- 
tion of the temperature at the critical density, which gives 
the background Ab = 2.3 far above the critical point. In 
Fig. 8 the data of the singular part AA = A — Ab are plot- 
ted as a functions of £. For £ < 10 our numerical data 
nicely agree with the theoretical linear response result 
Eq.19, which will be explained below. For £ > 10 the 
finite-size effect and the nonlinear response effect should 
be responsible for the saturation of the calculated A. 

The mode-coupling theory in Appendix B predicts the 
following behavior, 



A = \ B + -Lc p \n(L/0 
Attt] 



(19) 



where tj is the shear viscosity j25j. See also Appendix C 
for the linear response theory for heat flow. The singular 
part of the thermal conductivity is simply given by DtC p 
with Dt being the thermal diffusion constant 0. In 
terms of the isothermal compressibility Kt the coefficient 
A\ is written as 



( dp 



-7/4 



(20) 



To estimate A\ from the above expression, we calculated 
the shear viscosity 77 at T — 0.50 and n — 0.37 by two 



methods and obtained almost the same results. That is, 
(i) the time integral of the stress time-correlation func- 
tion 0,0] in the range < t < 100 gave 77 = 0.35 and (ii) 
the long time tail of the velocity correlation function gave 
7] = 0.33 (see Appendix D). If we use the latter result to- 
gether with Eqs.7 and 8, we are led to A\ = 0.035. In 
Fig. 8 the theoretical curve represents the second term on 
the right hand side of Eq.19 with this A\. It excellently 
agrees with the data before the saturation of A. 



D. Momentum and heat flux distributions at fixed 
density under heat flow 

Some characteristic features of the density fluctuations 
can be seen in the one-body density distribution function 
^(p) = (^(^ceii — p))i where n ce ii is the density in an ap- 
propriately chosen cell in the fluid and the average is 
taken over the thermal fluctuations and over the cells in 
the system [13 III 13 El • It is the probability distribu- 
tion of the coarse-grained density. Furthermore, in the 
presence of heat flow, we are interested in distributions 
of the momentum and heat fluxes at fixed density. They 
can give the correlations between these fluxes and the 
density within the same cell even if the cluster motion 
driven by heat flow is very small. 

First we coarse-grain the system to calculate ^(p). 
The interior region (— 0.4L < z < 0AL and < y < L) is 
divided into 10 x 10 rectangular subsystems. Let Mkit) 
(k = 1, • • ■ , 100) be the particle number in the fc-th cell 
at time t. After the time averaging in steady states, we 
obtained the distribution of Mk (t) for integer M as 



, 100 

P ( M ) = I7j7jE<<W^)>> 



(21) 



fc=i 



where 5m, M' is the Kronecker delta, and X)m=o P(M) = 
1 by definition. For each given density p = M/V^ e ll we 
define 



*(p) = V cea P(V cellP ) 



1 100 



(22) 



fc=i 



where V cc \\ = 0.8L 2 /100 is the cell volume and the second 
line is the expression in the continuum limit with (t) = 
Mk{t) /V cc \\. By definition we obtain 

/>oo />oo 

/ dpV(p) = l, / dppV{p) = n in , (23) 
Jo Jo 

where r7 in is the average density in the interior and n ln = 
n in our case. The second moment becomes 

^00 1 100 

^ dp(p-n in )H(p) = — jy,(n k (t)-n in ) 2 ). (24) 

In equilibrium, or if the heterogeneity along the heat flow 
is neglected, the second moment behaves as t, 2 ^ 71 /V ce n for 
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£ less than the cell length but as V^ d ^ fl ^ d for larger £ 
due to the finite-size effect. 

Now we consider the coarse-grained momentum and 
heat fluxes at fixed density. We calculate the following 
steady-state averages, 

^ 100 

J ^ p) = T,( J ok(t)HP ~ »*(*))) (25) 



, 100 

Jq(p) = yqqT? — E( J ot - ^w))' ( 26 ) 

ccl1 fe=l 

where ^(i) and J^, 2 : (t) are the z component of the 
space integral of the momentum density and that of the 
heat flux, respectively, within the fc-th cell (see (A. 2) in 
the Appendix A for their definitions). If they are di- 
vided by the cell volume V^di, they become the cores- 
grained densities, respectively. For simplicity, we may 
write ^(p) = (S(p — n)), J P (p) = (J z S(p ~ n)) and 
Jq{p) — (J®8(p — n)) regarding the dynamic variables 
involved as the coarse-grained quantities. The normal- 
ized quantities J p {p)/^{p) and Jq(p)/^(p) may be in- 
terpreted as the coarse-grained conditional average of the 
momentum density and that of the heat flux, respectively, 
under the condition of fixed density at p. If integrated 
over p, we obtain 

/>oo 

/ d P J p (p) = 0, (27) 
Jo 

/ d P J Q (p) = -Q, (28) 
Jo 

where Q is the average heat flux defined by Eq.18 in the 
interior. In Appendix C we will examine the expected be- 
havior of these quantities using the linear response theory 
for VT |H. 

In Fig. 9 we show the three quantities, ^(p), J P (p), and 
Jq(p), obtained from 10 independent runs. The temper- 
ature at z = L is Tl = 0.65 in (a) (upper plate), XL = 0.5 
in (b) (middle plate), and Tl = 0.48 in (c) (lower plate), 
with AT = 0.005 or (dT/dz) = 0.43 x 10~ 4 . As can 
be seen in Fig. 7, the calculated thermal conductivity is 
A = 5.96 in (a), 5.66 in (b), and 2.63 in (c). Salient fea- 
tures are as follows. 

(i) The density distribution ^(p) has a rather sharp peak 
in (a), a broad (still single) peak in (b), and double peaks 
in (c). We also calculated ^(p) in equilibrium at the same 
temperatures, which exhibits double flattened peaks for 
T = 0.5 and sharper double peaks for T = 0.48, so 
the double peak behavior emerges more conspicuously 
in equilibrium. Furthermore, as a complicating factor in 
heat flow, Fig. 5 indicates that the average density profile 
is considerably dependent on z in (b) and (c). 

(ii) The momentum distribution J p (p) is positive for p > 
0.37 and negative for for p < 0.37. This is consistent with 



the anti-symmetric behavior, J p (p) ~ Q(p — 0.37)^( j o), 
close to the criticality in Eq.(C.5) of Appendix C. Ev- 
idently, the liquid-like clusters move toward the higher 
temperature boundary, while the particles in the gas-like 
regions move toward the lower temperature boundary. 
However, notice that the high-density maximum is con- 
siderably sharper than the low-density minimum, which 
should arise from the gas-liquid asymmetry of the fluctu- 
ations |2^. In particular, for the case (b), the momentum 
density of the liquid-like regions is of order 10 -3 and the 
velocity is of order 3x 10~ 3 (in units of a /to = {e/m) 1 / 2 ). 
In this case we have £ ~ 18 and Dt ~ 0.1 so that 
the distance of the cluster motion within the life time 
i 2 /Dt ~3x 10 3 is estimated to be of order 10. 
(iii) The heat flux distribution function Jq(p) still ex- 
hibits considerable irregular behavior, but its negativity 
at any p is clear. Let us smooth out the curves; then, 
Jq(p) has a single minimum in (a) and double min- 
ima in (b) and (c). Thus, as T — > T c , heat is largely 
transported by the counterflow of the liquid-like clus- 
ters and the gas-like regions. Particularly in (c), the 
contribution from p = 0.37 becomes very small and 
the curve can be fairly fitted to the symmetric relation 

Jq(p) Q{p~0.37) 2 ^(p) in accord with Eq.(C. 6). The 

gas-liquid asymmetry is more suppressed for Jq(p) than 
for J p (p). 



IV. CONCLUDING REMARKS 

MD simulations have been performed on LJ near- 
critical fluids in two dimensions. In equilibrium the crit- 
ical properties obtained are presented in Figs. 1-3. The 
main results under heat flow are summarized as follows. 

(i) We have calculated the average density and temper- 
ature profiles in a steady state in Fig. 4, where they are 
fairly fitted to linear lines and satisfy Eq.16. The density 
deviation is much enhanced than that of the temperature 
and the average pressure remains homogeneous. 

(ii) We have obtained critical enhancement of the ther- 
mal conductivity for various T close to T c in Figs. 7 and 
8 in good agreement with the mode-coupling prediction 
in Eq.19 derived in Appendix B. 

(iii) We have calculated the one-body density distribu- 
tion ^(p), the momentum distribution J p (p), and the 
heat flux distribution Jq(p) defined by Eqs.21, 24, and 
25. Fig. 9 demonstrates the cluster convection mecha- 
nism, which is briefly summarized in the introduction 
and supported in Appendix C in the linear regime. 

(iv) The cluster convection is a natural consequence of 
the irreversibility in heat conduction, while the density 
increase near the cooler boundary in Fig. 6 arises from the 
simple thermodynamics under homogeneous pressure in 
Eq.16. These two effects are not contradictory with each 
other in view of the fact that the distance of cluster con- 
vection is very short. 

The following problems could be mentioned as future 
subjects of nonequilibrium MD simulations. 
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(i) When the boundary wall is heated with a fixed cell vol- 
ume, sound waves emitted from the boundary can cause 
rapid adiabatic heating throughout the cell (the piston ef- 
fect) |27l l28l| . We should examine how this phenomenon 
starts in the early stage on the acoustic time scale |2( . 

(ii) Heat conduction in two-phase near-critical fluids be- 
low T c has been little examined in the literature 0. For 
example, we should examine how a gas-liquid interface re- 
acts to applied heat flow, where latent heat transport can 
be crucial in the presence of convection. Interestingly, gas 
bubbles in liquid migrate toward the warmer boundary 
in heat flow owing to the Marangoni effect p9j.@ 
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Appendix A: Microscopic Expressions 

We introduce the momentum density, 



The Green-Kubo formula for the thermal conductivity 
reads 



A = 



The J^it) in Eq.18 is the z component of the total heat 
flux in the interior, 

J?(t) = f drJ Q (r,t). (A.5) 

J interior 

In Eqs.24 and 25 the space integrals are within small 
subsystems. 

Appendix B: Mode-Coupling Theory 

In the critical dynamics of simple fluids the gross vari- 
ables include the long wavelength parts (with wave num- 
bers in the region 9 « u -1 ) of the energy density e, the 
particle density n, and the momentum density J. The 
heat flux density J (r, t) in (A. 3) has been approximated 
as a sum of a product of the gross variables and a random 
part in the form 0,11110, 

J Q (r, t) = —Ss(r, t)J(r, t) + J%(r, t). (B.l) 

m 

The Ss is the fluctuating entropy deviation (per particle) 
defined by 



1 



fc R T 2 



dt / dr ( Jf (r,t)J? (0,0)}. (A.4) 



J(r,t) =Y^rnVi{t)8{r - r 4 (t)) 



(A.l) 



and the energy current density J c (r,t). The micro- 
scopic expression for the latter quantity is rather com- 
plicated Q. Let us consider its space integral Jo(t) — 
f v drJ c (r, t) in a subsystem with volume V\ containing 
many particles. It may be approximated as 



o£'£^( r «) — ( v « ■ r v) r V> ( A - 2 ) 



where r 4 = r 4 (t) and Vi — Vi(t) are the position and ve- 
locity of the i-th particle (the time t being suppressed in 
(A.2)), Tij = Vi — Tj, 4>'(r) = d<p(r)/dr, and the summa- 
tion ^T. is over the particles contained in the subsystem 
under consideration. Here the pair interactions between 
the particles inside and outside the subsystem are not 
precisely accounted for. 

The microscopic heat flux density is defined by 

J Q {r,t) = J c (r,t) - [(e + p)/n]J(r,t), (A.3) 

where e, p, and n are the average energy, pressure, and 
density, respectively. This current satisfies the orthogo- 
nal property J dr(J®{r,t) ■ J(r', i)) = in equilibrium. 



8s(r,t) 



1 

nT 



8e(r,t) 



P 



8n(r, t) 



(B.2) 



in terms of the deviations of the energy density e and the 
number density h The e can be defined microscopically 
using the particle positions and velocities 0, |24(- The 
first term on the right hand side of (B.l) evolves slowly 
in time and gives rise to the singular part of the thermal 
conductivity AA when substituted into (A.4). In 2D the 
mode-coupling calculation yields the following integral 
over the wave vector q, 



AA = 



k B T 
2rj 



dq 1 

(2^)2^2 



C P (q), 



(B.3) 



where r\ is the shear viscosity |25j and C p (q) = 
fcg 1 n 2 (|sqi| 2 ) is the variance of the entropy fluctuation 
with sq being the Fourier component. See Appendix C 
for another derivation of AA from the linear response. As 
far as the most singular part is concerned, we may set 



This yields 



8s = -n- 2 (dp/dT) c ^8n. 



C p (q) S (dp/dT) 2 cx S(q)/k B n, 



(B.4) 



(B.5) 



in terms of the structure factor S(q) in Eq.5 2]. The 
long wavelength limit C p — lim^o C p (q) is the usual iso- 
baric specific heat per unit volume behaving as in Eq.10. 
Note that the integral (B.3) is logarithmically divergent 



7 



at small q, so we obtain the expression Eq.19. On the 
other hand, the second term on the right hand side of 
(B.l) relaxes rapidly and gives rise to the background 
thermal conductivity Ab- 

Appendix C: Linear Response to Temperature 
Gradient near the Gas-Liquid Critical Point 

Here we consider the linear response theory with re- 
spect to a temperature gradient a = VT (along the z 
axis) in a steady heat-conducting state in the absence 
of macroscopic velocity field 26] . To pick up the singu- 
lar contribution near the gas-liquid critical point we may 
approximate the heat flux by (T/m)Ss(r,t)J(r,t) from 
(B.l). Then the linear response of any dynamic variable 
B(r, t) to a can be written as 



x (S(p — h(r))Sn(r')), 



(C.3) 



x{6(p-h(r))8n(r)Sn(r')). (C.4) 
We notice that these quantities depend on the cell volume 
V^cii- If the cell length £ cc \\ = Vllf is shorter than the 
correlation length £, we estimate J p (p) as 



J p(p) ~ — 



6(B) 



—a 



where n- ln is the average density. If £ ce n is longer than £, 
we divide the cell into subsystems with length £ and find 
that J p (p) is given by (C.5) with ^ ell being replaced by 



mk B T J 



poc t> y v ' ° — •J \ / ecu o <j 

/ dt dr'(B(r,t)5s{r',0)J(r\0)) £ 2 - Next notice that the integral /dp J Q (p) is equal to 
7n J —(A\)dT/dz from (C.4) which is in accord with Eq.28 



dt dr'(B(r,0)Ss(r',t)J(r',t)). (C.l) 



-(A\)dT/dz from (C.4) which is in accord with Eq.28 
for AA = A. Accounting for this sum rule we thus expect 



rnksT J 



From the first to second line use has been made of the 
time-reversal relation (A(t)B(0)} = (B(t)A(G)) where A 
and B are the time-reversed variables. For example, 
J = J. Furthermore, on the second line, we may re- 
place 5s(r', t) by 5s(r', 0) because the relaxation time of 
J(r',t) due to the shear viscosity -q is much faster than 
that of Ss(r',t). Then the time integral may be per- 
formed to give 



Jq(p) = ~A Q {p-n m )H{p) 



dT 

dz ' 



(C.6) 



8(B) 



x(S(r)S S (r')J,(r")>, 



(C.2) 



where the equal-time correlation is involved and the time 
dependence is hence suppressed. The %j(r) is the Oseen 
tensor whose Fourier transformation is %j (q) = (5ij — 
QiQj /q 2 )/il < l 2 ■ I n 3D it follows the well-known expression 
Tij(r) = (Sij + XiXj/r 2 )/8TTr]r. 

In (C.2), if we set B = (T/m)SsJ z and use the equilib- 
rium relation (Ji(r)Jj(r')) = UrT p5ij8(r — r') : we repro- 
duce the mode-coupling expression for the singular part 
of the thermal conductivity (given by (B.3) in 2D) in the 
form 8(B) = -AXdT/dz. Next let us set B = J z S(p - ti) 
and jQS(p-n) where the dynamic variables J z , jQ, and 
n are the coarse-grained quantities averaged in appropri- 
ate cells. Then we obtain J p (p) and Jq(p) in Eqs.24 and 
25 expressed as 



for £ cc \\ <C £. The coefficient Aq is determined from the 
normalization condition Eq.28. The estimations (C.5) 
and (C.6) are consistent with the data in Fig. 9. 

In addition, Eq.l in the introduction follows if we as- 
sume V£ ~ J p (p)/mn 1 $(p) in (C.5) by setting £ cc \\ ~ £ 
and p — n m ~ £ with the aid of the exponent relation 
2/3 = (d-2 + f])v 0. Note that J p {p)/^{p) represents 
the average momentum density at density p. 

Appendix D: Diffusion in Two Dimensions 

In two dimensions the flux-time correlation functions 
for the transport coefficients have a long time tail re- 
laxing as 1/t, giving rise to a logarithmic singularity (if 
integrated over time) [25|. The simplest example is the 
diffusion constant D of a tagged particle. It is the time 
integral of the velocity time-correlation function, 



G(t) 



1 N 

— ^2(Vi(to+t)'Vi{to)). 



(D.l) 



dp 
dT 



The long time tail of G(t) is theoretically given by 
{k^T /Smf)/t if the kinetic viscosity rj/mn is much larger 
than D. By taking the average over to in a time interval 
of 5 x 10 4 , we obtained ^dt'G{t') = 0.17 + 0.059 \ogt 
for t > 1, leading to k B T/8-KT] = 0.059. Note that the 
kinetic viscosity is close 1 and is considerably larger than 
the diffusion constant in our system. 



[1] J.V. Sengers and J.M.H. Levelt Sengers, in Progress in [2] A. Onuki, Phase Transition Dynamics (Cambridge Uni- 
Liquid Physics, edited by C.A. Croxton (Wiley, Chich- versity Press, Cambridge, 2002). 

ester, England, 1978). [3] L.P. Kadanoff and J. Swift, Phys. Rev. 166, 89 (1968). 



8 



[4] K. Kawasaki, in Phase Transition and Critical Phenom- 
ena, edited by C. Domb and M.S. Green (Academic, New 
York, 1976), Vol.5A; Ann. Phys. (N.Y.) 61, 1 (1970). 

[5] E.D. Siggia, B.I. Halperin, and P.C. Hohenberg, Phys. 
Rev. B 13, 2110 (1976). 

[6] J. P. Hansen and I.R. McDonald, Theory of Simple Liq- 
uids, (Academic Press, London, 1986). 

[7] M.P. Allen and D.J. Tildesley, Computer Simulation of 
Liquids (Clarendon Press, Oxford, 1987). 

[8] R. Vogelsang, C. Hoheisel, G. V. Paolini, and G. Ciccotti, 
Phys. Rev. A 36, 3964 (1987); R. Vogelsang, C. Hoheisel, 
and G. Ciccotti, J. Chem. Phys. 86, 6371 (1987). 

[9] P. J. Gardner, D. M. Heyes, and S. R. Preston, Mol. 
Phys. 73, 141 (1991). 
[10] G. V. Paolini, G. Ciccotti, and C. Massobrio, Phys. Rev. 

A 34, 1355 (1986). 
[11] D. J. Evans and S. Murad, Mol. Phys. 68, 1219 (1989); 
B. Y. Wang, P. T. Cummings, and D. J. Evans, Mol. 
Phys. 75, 1345 (1992). 
[12] B. Hafskjold, T. Ikeshoji, and S. K. Ratkje, Mol. Phys. 
80, 1389 (1993); T. Ikeshoji and B. Hafskjold. ibid. 81, 
251 (1994). 

[13] T. Ohara, J. Chem. Phys. Ill, 9667 (1999). 
[14] T. Ohara, J. Chem. Phys. Ill, 6492 (1999). 
[15] J. A. Barker, D. Henderson, and F.F. Abraham, Physica 

106A, 226 (1981). 
[16] R.R. Singh, K.S. Pitzer, J.J. de Pablo, and J.M. Praus- 

nitz, J. Chem. Phys. 92, 5463 (1990). 
[17] B. Smit and D. Frenkel, J. Chem. Phys. 94, 5663 (1990). 
[18] S. Jiang, and K.E. Gubbins, Molec. Phys. 86, 599 (1995). 
[19] M. Rovere, J. Phys.; Condens. Matt. 5, B193 (1993); 

M. Rovere, P.Nielaba, and K. Binder, Z. Phys. 90, 215 

(1993). 

[20] H. Luo, G. Ciccotti, M. Mareschal, M. Meyer, and B. 
Zappoli, Phys. Rev. E 51, 2013 (1995). 

[21] S. Nose, Molec. Phys. 52, 255 (1983). 

[22] On the critical isochore {dp/dT)n/(dp/dT) cx - 1 is of 
order C v /C p ~ (T/T c - 1) 7 " Q @, where 7 = 7/4 and 
q = in two dimensions. 

[23] In Ising spin systems we may equally define clusters 
with up-spins and those of down-spins. In our simula- 
tions, however, only the denser clusters (n > n c ) are well 
defined as in Fig. 5. The symmetry between the higher- 
density and lower-density critical fluctuations should be 
attained much closer to the criticality. 

[24] N.B. Wilding, J. Phys. Condens. Matt. 9, 585 (1997). 

[25] B. J. Alder and T. E. Wainwright Phys. Rev. A 1, 18 
(1970); J. R. Dorfman and E. G. D. Cohen, Phys. Rev. 
Lett. 25, 1257 (1970); T. E. Wainwright, B. J. Alder, 
and D. M. Cass, Phys. Rev. A 4, 233 (1971). The shear 
viscosity also has a small logarithmic singularity r\ = 770 + 
^4.,, log(L/o-) with A v /r/o ~ mnk-BT /9,-KriQ being small. 

[26] I. Procaccia, D. Ronis, M.A. Collins, J. Ross, and I. Op- 



penheim, Phys. Rev. A 19, 1290 (1979). 
[27] A. Onuki, H. Hao and R. A. Ferrell, Phys. Rev. A 41, 
2256 (1990). 

[28] Y. Garrabos, M. Bonetti, D. Beysens, F. Perrot, T. 
Frohlich, P. Carles, and B. Zappoli, Phys. Rev. E 57, 
5665 (1998). 

[29] D. Beysens, Y. Garrabos, V. S. Nikolayev, C. Lecoutre- 
Chabot, J. -P. Delville, and J. Hegseth, Europhys. Lett., 
59, 245 (2002). 



100 



10 



(JO 



0.1 




0.05(~27r/L) 



_l I 1 ' I I I I 



0.01 



0.1 



10 



FIG. 1: The structure factor S(q) at@ n = 0.37 for T = 
0.65 (short-dashed line at bottom), 0.51, 0.5, 0.495, and 0.49 
(dashed line on top). A line with a slope of —7/4 is included 
as a guide. 
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FIG. 2: The correlation length £ vs the density at various 
temperatures obtained from the structure factor. 
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FIG. 3: The isothermal compressibility Kt vs the density at 
various temperatures obtained from the structure factor. 
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FIG. 4: Simulation cell under heat flow composed of cooling, 
heating, and interior regions |T^I . 
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FIG. 5: Snapshot of the particle configuration at n = 0.37 in 
a steady state with Tl = 0.50 and Th = 0.505. The horizontal 
bars at the vertical box lines mark the boundary between the 
interior region and cooling or heating region. 




FIG. 6: Steady-state temperature and density profiles in the 
z-direction obtained by the time average. The solid line has a 
slope of AT/(9L/10) with AT = 0.005, while the dotted line 
has a slope of -na p AT/(9L/10) in Eq.16. 
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FIG. 7: The thermal conductivity A calculated from Eqs.16 
and 17 at n = 0.37 for T = 0.7, 0.65, 0.6, 0.52, 0.51, 0.5, 
0.495, and 0.49. The bold dashed line is a view guide. The 
width of each error bar is twice of the variance of 10 data 
values corresponding to 10 independent runs. 




FIG. 8: The singular part of the thermal conductivity AA as 
a function of £ on logarithmic scales. The solid line is the 
second term in Eq.19 with A\ — 0.035. 
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FIG. 9: (a) Density distribution function \&(p) (right scale), 
momentum distribution J p (p), and heat-flux distribution 
Jq(p) (left scale) obtained at n = 0.37 for TL = 0.65 in (a), 
T L = 0.50 in (b), and T L = 0.48 in (c). 



